//Difinition of interval:
typedef float2 interval;

interval iadd(interval a, interval b) {
	//return interval( add(a.x, b.x), add(a.y, b.y) );
	return interval( (a.x+b.x), (a.y+b.y) );
}
interval imul(interval a, interval b) {
	float4 lh = a.xxyy * b.xyxy;
	return interval(min(lh.x, min(lh.y, min(lh.z, lh.w))),
					max(lh.x, max(lh.y, max(lh.z, lh.w))));
}
interval ircp(const float inf, interval i) {
	return ( (i.x <= 0 && i.y >= 0) ? interval(-inf, inf) : 1/i.yx );
}

//Difinition of raf:
typedef float3 raf;
raf interval_to_raf(interval i){
	raf r;
	r.x = (i.y + i.x);
	r.y = (i.y - i.x);
	r.xy *= .5; r.z = 0;
	return r;
}

float raf_radius(raf a){
	return abs(a.y) + a.z;
}

interval raf_to_interval(raf a){
	const float rad = raf_radius(a);
	return interval(a.x - rad, a.x + rad);
}

raf raf_add(raf a, raf b){
	return (a + b);
}

raf raf_mul(raf a, raf b){
	raf r;
	r.x = a.x * b.x;
	r.y = a.x*b.y + b.x*a.y;
	r.z = abs(a.x*b.z) + abs(b.x*a.z) +
		  raf_radius(a)*raf_radius(b);
	return r;
}

bool raf_contains(raf a, float t){
	return (abs(t-a.x)<raf_radius(a));
}

raf raf_add_const(raf a, float t){
	raf re=a;
	re.x = re.x + t;
	return re;
}

raf evaluate_raf(raf X, raf Y, raf Z/*, float w, float nan, float inf*/){
	//raf re=raf_mul(X,X);	  return raf_add(raf_mul(re, raf_add_const(re, -1)), raf_add(raf_mul(Y,Y), raf_mul(Z,Z)));	//case 17 Wrong
	raf re = raf_mul(X,X);  return raf_add(raf_add(raf_mul(re, raf_add_const(re, -1)), raf_mul(Y,Y)), raf_mul(Z,Z));  //case 17 Correct
	//return raf_add_const(raf_add(raf_add(raf_mul(X,X),raf_mul(Y,Y)), raf_mul(Z,Z)), -1);								//Sphere
}


float eps;
int width;
int height;
#define _FX_ZERO_EPS eps
float traverse(float3 penter, float3 pexit){
			   //, float w, float max_depth, float nan, float inf

	const float3 org = penter;
	const float3 dir = pexit-penter;
	interval t = interval(0,1);
	raf F, ix, iy, iz;

	//rejection test
	//it = interval_to_raf(t);
	//ix = raf_add(org.x, raf_mul(it, dir.x));
	//iy = raf_add(org.y, raf_mul(it, dir.y));
	//iz = raf_add(org.z, raf_mul(it, dir.z));
	//F = evaluate_raf(ix, iy, iz, w, nan, inf);
	ix = interval_to_raf(interval(org.x+dir.x*t.x, org.x+dir.x*t.y));
	iy = interval_to_raf(interval(org.y+dir.y*t.x, org.y+dir.y*t.y));
	iz = interval_to_raf(interval(org.z+dir.z*t.x, org.z+dir.z*t.y));
	F = evaluate_raf(ix, iy, iz);
	
	if (raf_contains(F, 0)){
		int d=0;
		float tincr = .5;
		const int dlast = log2(length(dir)/eps);
		//main loop
		for(;;){
			t.y = t.x + tincr;
			
			ix = interval_to_raf(interval(org.x+dir.x*t.x, org.x+dir.x*t.y));
			iy = interval_to_raf(interval(org.y+dir.y*t.x, org.y+dir.y*t.y));
			iz = interval_to_raf(interval(org.z+dir.z*t.x, org.z+dir.z*t.y));
			F = evaluate_raf(ix, iy, iz);

			if (raf_contains(F, 0)){
				if (d==dlast){ return t.x; /*hit*/}
				else{ tincr *= .5; d++; continue; }
			}
			else{
				t.x = t.y;

				float fp = frac(.5*t.x/tincr);
				if (fp < 1e-8){
				//back-recursion
					for(int j=0; j<=dlast; j++){
						tincr *= 2;
						d--;
						fp = frac(.5*t.x/tincr);
						if (d==-1 || fp > 1e-8) break;
					}
					if (d==-1) break;//No hit
				}
				//else: just step forward with step=tincr.
			}
		}
	}
	return -1; //no hit
}

void RayCastingPS(float2 vPos:VPOS,
				out float4 color:COLOR0)
{
	float4 dir = float4(2*vPos.x/width-1.0f, 1.0f-2*vPos.y/height, 1.0f, 1.0f);
	dir = mul(dir, g_mInvWVP);
	dir = dir/dir.w;
	dir.xyz = dir.xyz-eye;
	
	float tmin=-3.402823466e+38F;
	float tmax=3.402823466e+38F;
	if(dir.x>_FX_ZERO_EPS){
		tmin = max(tmin, (box0.x   -eye.x)/dir.x);
		tmax = min(tmax, (boxLast.x-eye.x)/dir.x);
	}
	else if(dir.x<-_FX_ZERO_EPS){
		tmin = max(tmin, (boxLast.x-eye.x)/dir.x);
		tmax = min(tmax, (box0.x   -eye.x)/dir.x);
	}

	if(dir.y>_FX_ZERO_EPS){
		tmin = max(tmin, (box0.y   -eye.y)/dir.y);
		tmax = min(tmax, (boxLast.y-eye.y)/dir.y);
	}
	else if(dir.y<-_FX_ZERO_EPS){
		tmin = max(tmin, (boxLast.y-eye.y)/dir.y);
		tmax = min(tmax, (box0.y   -eye.y)/dir.y);
	}

	if(dir.z>_FX_ZERO_EPS){
		tmin = max(tmin, (box0.z   -eye.z)/dir.z);
		tmax = min(tmax, (boxLast.z-eye.z)/dir.z);
	}
	else if(dir.z<-_FX_ZERO_EPS){
		tmin = max(tmin, (boxLast.z-eye.z)/dir.z);
		tmax = min(tmax, (box0.z   -eye.z)/dir.z);
	}
	if(tmin>tmax||tmax<0)	clip(-1);


	float hit=traverse(eye+tmin*dir.xyz, eye+tmax*dir.xyz);
	if(hit<0)	clip(-1);


	float3 Pos = eye+(tmin+hit*(tmax-tmin))*dir.xyz;
	color.rgb=Pos; color.a = 1.0; return;

	float3 dx = ddx(Pos);
	float3 dy = ddy(Pos);
	float3 Norm = cross(dx, dy);
	Norm = normalize(Norm);

	float3 matDiff = ChooseMaterial(Pos, Norm);
	color.rgb = shading(Pos, Norm, matDiff, nLight);
	color.a = 1.0f;
}

technique RayCasting
{
	pass p0
	{
		VertexShader = null;
		PixelShader  = compile ps_3_0 RayCastingPS();
	}
}